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Abstract 

In  many  spray  applications,  it  is  important  to 
know  the  size  and  velocity  distribution  of  the  drops. 
Conventional  particle  tracking  techniques  can  require 
prohibitively  large  computational  times,  especially  in 
regions  of  low  droplet  number  density  or  when  detailed 
statistics  are  desired.  In  this  paper,  we  develop  a 
method  to  compute  the  statistics  directly  by  solving  a 
series  of  moment  transport  equations.  A  maximum 
entropy  model  is  used  to  close  higher-order  moments 
appearing  in  the  equations.  Solution  of  these  equations 
gives  not  only  the  transported  moments  of  the  spray, 
but  also  the  maximum  entropy  probability  distribution 
function  from  which  further  statistics  can  be  obtained. 
The  method  has  been  tested  on  a  quasi-one-dimensional 
spray  problem  to  assess  its  viability.  Submodels  which 
account  for  the  effects  of  the  gas  on  the  drops,  including 
turbulence  modification  and  correlation  between  the  gas 
and  drop  velocities,  are  incorporated.  Results  for 
expected  quantities  are  in  good  agreement  with  the 
solution  from  a  particle  tracking  simulation. 

Introduction 

Understanding  and  modeling  multiphase  flows  have 
become  of  primary  interest  in  a  wide  variety  of 
applications.  Of  particular  importance  are  spray  flows 
such  as  those  associated  with  liquid  fuel  injectors, 
industrial  coating  processes,  and  agricultural  sprays.  A 
spray  flow  can  be  defined  as  that  regime  downstream  of 
an  injector  where  a  liquid  column  or  sheet  has  broken 
up  and  atomized,  but  the  resultant  drops  continue  to 
have  some  mean  velocity  relative  to  the  gas  phase,  as 
opposed  to  an  aerosol  where  there  is  no  mean  slip 
velocity.  Both  the  liquid  phase  and  the  gas  phase 
continue  to  interact  dynamically,  exchanging  mass. 


momentum,  and  energy.  This  evolution  can  have  a 
significant  impact,  for  example,  on  the  combustion 
process  inside  a  liquid  rocket  engine  or  on  the  coating 
produced  by  a  spray  application  system.  For  this 
reason,  it  is  necessary  to  have  a  full  understanding  of 
the  physics  of  spray  flows  and  to  be  able  to  accurately 
predict  not  only  the  distribution  of  the  drops  and  their 
behavior,  but  also  the  dynamics  of  the  gas  phase,  that 
is,  the  characteristics  of  the  combined  two-phase  flow. 

Probably  the  most-used  approach  in  studying  spray 
flows  has  been  the  Lagrangian-Eulerian,  or  particle 
tracking  method.1,2,3,4  In  this  approach,  the  gas  phase  is 
predicted  by  solving  the  time-dependent,  Reynolds- 
Averaged  Navier-Stokes  (RANS)  equations  with  a 
suitable  turbulence  model  and  appropriate  exchange 
terms.  Drops  are  stochastically  injected  into  the  gas  and 
the  drop  trajectories  are  computed  by  integrating  a 
Lagrangian  equation  of  motion. 

While  the  Lagrangian-Eulerian  approach  has 
provided  useful  information  in  many  applications,  it  has 
some  significant  drawbacks.  For  instance,  data  arising 
from  any  simulation  must  be  post-processed.  If  the 
quantity  of  interest  is  the  mean  number  density  in  each 
grid  cell,  this  may  not  pose  a  problem.  However,  what 
if  we  are  interested  in  the  mean  number  density  of  drops 
in  a  particular  size  class?  The  averaging  procedure 
might  be  simple,  but  the  required  computational  time  is 
increased  because  a  sufficient  number  of  drops  must 
pass  through  the  cell  to  provide  a  data  set  large  enough 
for  a  meaningful  average.  As  the  quantity  of  interest 
becomes  more  specific,  the  necessary  computation  time 
becomes  more  prohibitive. 

An  alternative  approach  which  does  not  involve 
simulation  is  to  compute  the  evolution  of  a  probability 
density  function  (PDF)  describing  the  drops.  Williams5 
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was  the  first  to  derive  a  transport  equation  for  a  drop 
PDF,  called  the  spray  equation,  analogous  to 
Boltzmann’s  equation  for  molecules.  However,  because 
obtaining  a  solution  to  the  spray  equation  is  quite 
difficult,  Williams  solved  a  simplified  version  of  the 
equation  for  a  steady,  one-dimensional  spray  flow  where 
the  gas  properties  were  known  in  advance.  Other 
solutions6,7,8,9,10  have  been  attempted,  but  with  limited 
success.  Because  of  the  high-dimensionality  of  the 
spray  equation,  computations  could  only  be  done  with  a 
coarse  discretization  of  the  phase  space. 

A  different  approach  to  solving  the  spray  equation 
was  taken  by  Tambour.11  He  only  considered  the  size 
distribution,  but  instead  of  transporting  the  PDF,  he 
discretized  the  size  axis  and  derived  a  set  of  sectional 
equations  within  each  size  bin.  The  advantage  of  this 
method  is  that  it  reduces  the  dimensionality  of  the 
equation  to  only  physical  space,  though  a  set  of 
transport  equations  for  each  bin  has  to  be  solved. 

The  purpose  of  this  paper  is  to  extend  the  PDF 
approach  and  obtain  a  complete  description  of  a  spray 
flow  by  computing  the  evolution  of  its  PDF  along  with 
the  gas  flow  in  which  it  is  embedded.  This  will  be 
accomplished  by  deriving  and  solving  a  set  of  moment 
transport  equations  for  average  quantities  of  interest, 
such  as  the  mean  drop  velocity  and  diameter,  and 
developing  an  appropriate  closure  model. 

Governing  Equations 

In  deriving  the  equations  that  will  describe  the 
evolution  of  our  system,  we  will  make  use  of  ensemble 
averaging.  This  is  an  averaging  procedure  that  is 
independent  of  any  time  or  length  scales.  We  will  also 
restrict  ourselves  to  single-point  statistics.12  Before 
constructing  a  probability  density  function  from  which 
average  quantities  can  be  obtained,  it  is  necessary  to 
discuss  some  assumptions  about  the  spray.  A  single 
drop  can  be  described  by  any  number  characteristics,  or 
marks,  such  as  its  size,  velocity  or  temperature.  The 
trade-off,  though,  is  in  increasing  complexity  as  more 
marks  are  used.  Therefore,  it  is  often  sufficient  to 
choose  a  relatively  small  set  of  characteristics  that 
contains  the  information  of  interest.  This  set  of 
characteristics  prescribes  a  state  vector  a  whose 
elements  define  the  axes  of  a  hyperspace,  called  a- 
space,  through  which  the  drop  can  move.  For  example, 
if  the  chosen  characteristics  of  the  drop  are  its  position 
in  physical  space  and  its  diameter,  then 
a  ~  [x^x2,x3,(p)  and  a -space  is  a  four-dimensional 
hyperspace  whose  axes  are  the  three  position  coordinates 
jcpjc2,jc3  and  the  diameter  (j).  As  the  drop  moves 
through  physical  space  or  as  it  evaporates  or  condenses, 
its  location  in  a  -space  changes. 


Because  we  are  dealing  with  point  particles,  volume 
displacement  effects  are  negligible,  and  we  are  able  to 
treat  any  interphase  exchange  as  if  the  drop  were  alone 
in  a  locally  uniform  gas  field. 

Under  these  assumptions,  a  probabilistic 
description  of  the  spray  can  be  developed,12  defining 
as  the  probability  density  of  finding  a  single 
drop  at  a  point  (3c,  t).  It  is  also  the  expected  number 
density  of  drops  at  x  and  t.  This  dual  interpretation  is 
important  since  we  are  generally  more  interested  in  the 
expected  number  density  rather  than  the  probability 
density.  The  probabilistic  description  also  defines  a 
function  f(a'\  x,t)f  the  probability  density  of  finding  a 
particular  diameter  and  velocity,  conditioned  on  a  drop 
existing  at  x  and  t ,  where  a V  -space  is  the  subset  of  a  - 
space  excluding  the  spatial  coordinate.  By  combining 
these  two  functions,  we  define  the  unnormalized  single 
particle  probability  density  function  for  the  spray 

F{a\t)da  =  X[x\t)  f{a' \x,t)da'dV  (1) 

F{a\t)  is  unnormalized  because  A(x;f),  also  an 
unnormalized  PDF,  determines  how  much  spray  is 
present  at  x  and  t. 

If  a  =  (3c,0,v,T),  a  transport  equation  for  F[a;t) 
can  now  be  derived.13,14,15 

¥+V1-(Fv)+Vv-(F5) 

3  (2) 

+^(FO)  +  ^(F0)  =  ^+il+^° 

Equation  (2)  is  known  as  the  spray  equation.  It 
describes  the  evolution  of  the  probability  distribution 
function  F(jc,0,v,7;f)  through  joint  physical, 
diameter,  velocity,  and  temperature  space  and  includes 
source  terms  accounting  for  binary  collision  S2 ,  unary 
breakup  5, ,  and  zero-body  events  50  such  as  nucleation 
and  complete  vaporization.  The  spray  equation  has  not 
generally  been  viewed  as  a  practical  way  of  predicting 
spray  flows.  It  is  not  an  ordinary  evolution  equation  in 
the  sense  that  the  quantity  of  interest  is  only  being 
transported  through  physical  space,  but  also  through 
diameter,  velocity,  and  temperature  space.  A  numerical 
solution  requires  a  full  discretization  of  this  hyperspace. 
This  makes  obtaining  a  direct  solution  difficult. 

Despite  the  difficulties,  attempts  have  been  made  at 
solving  the  spray  equation.5,6,16,17  Some  of  these 


*  For  clarity  when  denoting  a  PDF,  those  arguments  over  which  the 
function  is  a  density  will  be  listed  first,  followed  by  a  semi-colon, 
followed  by  arguments  that  are  parameters  of  the  function. 
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solutions  make  use  of  an  assumed  form  for  the  PDF  F, 
which  leads  to  a  set  of  transport  equations  for  the 
parameters  that  describe  the  assumed  form.  The  primary 
problem  with  this  approach  is  that  there  is  no  guarantee 
that  the  spray  will  conform  to  the  assumed  distribution 
throughout  its  evolution,  or  even  that  F  can  always  be 
described  analytically. 

However,  even  the  most  general  PDF  can  be 
described  by  an  infinite  set  of  moments.  For  instance,  a 
one-dimensional  distribution’s  mean  determines  its 
location  on  the  axis,  its  variance  is  a  measure  of  its 
width,  the  third  central  moment  determines  the  amount 
of  skewness,  etc.  The  value  of  each  moment  affects  the 
shape  of  the  distribution  in  some  way,  but  notice  that 
each  higher-order  moment  affects  the  shape  in  a  less 
drastic  manner  than  the  next  lower-order  moment. 
Thus,  it  is  logical  to  assume  that  the  lowest-order 
moments  carry  sufficient  information  to  reasonably 
approximate  the  shape  of  the  distribution  as  a  whole.  If 
it  is  possible  to  obtain  a  function  for  this  shape  from 
these  low-order  moments,  we  can  then  approximate  all 
the  other  information  included  in  the  actual  PDF. 
Then,  by  solving  a  set  of  transport  equations  for  these 
moments,  we  can  determine  how  the  approximate 
distribution  evolves  in  time. 

Splitting  the  Spray  Equation 

Before  continuing,  we  make  some  simplifying 
assumptions  about  the  spray.  First,  we  will  limit  our 
description  of  a  drop  to  its  spatial  location  x ,  its 
velocity  v,  and  its  diameter  0.  Since  we  are  not 
including  the  drop  temperature,  we  will  only  be 
considering  nonvaporizing  sprays.  The  number  density 
of  drops  will  be  sufficiently  small  such  that  droplet 
collisions  can  be  neglected.  Finally,  unary  breakup  and 
zero-body  source  terms  will  not  be  considered.  These 
assumptions  do  not  limit  the  validity  of  the 
model-what  is  left  out  can  be  included  at  any  point  later 
on.  We  make  them  here  for  the  sake  of  simplicity. 

Upon  applying  these  assumptions  to  the  spray 
equation  (2),  we  get 

r)  F 

— +V,-(Fv)+Vy.(F3)  =  0  (3) 

Substituting  equation  (1)  in  the  form 

F(0,  v,  x ;  t )  =  X(x ;  v;x,f)  (4) 

into  (3)  and  carrying  out  the  product  rule  gives 


A^  +  Z^  +  AV,  (/v)+ /v- V, A 
dt  dt  K  ’  (5) 

+AVv  •  (fa) + /5  •  Vv  A  =  0 

By  integrating  this  equation  over  all  velocity  and 
diameter  space,  i.e.  a' -space,  to  eliminate/,  we  obtain 

f^+V  -(A(v»  =  0  (6) 

where 

(v)  =  J vf(cc')da'  (7) 


is  the  ensemble-averaged,  or  expected,  drop  velocity. 
Equation  (6)  describes  the  evolution  of  the  expected 
number  density  of  drops  through  space  and  time. 
Notice  that  there  is  no  reference  to  a' -space  in  this 
equation.  The  integration  has  reduced  the  dimension  of 
the  equation  from  seven  (three  spatial,  one  diameter,  and 
three  velocity  coordinates)  to  three. 

Multiplying  equation  (6)  by  /  and  subtracting  the 
result  from  equation  (5),  we  get,  after  rearranging, 


f  +  Vj.(/v)  +  Vv.(/5)  = 
/[V,<v>+«v>-v).V>A] 


(8) 


This  equation  describes  the  evolution  of  the  conditional 
probability  density  function  through  a '  -space.  If  we 
interpret  (8)  as  a  flow  equation,  then  the  first  term 
accounts  for  the  accumulation  of  probability.  The  next 
two  terms  represent  the  advection  of  probability  through 
physical  and  velocity  space,  respectively.  Additionally, 
there  are  two  source  terms. 

The  first  term  on  the  right  side  represents  the 
production  of  probability  due  to  gradients  in  the  mean 
drop  velocity.  Suppose  that  there  is  a  positive  gradient 
in  the  mean  velocity  of  drops  moving  through  a  small 
volume  of  physical  space.  Because  of  the  gradient, 
drops  will  be  leaving  the  volume  at  a  greater  rate  than 
they  enter.  This  represents  a  loss  of  probability. 
However,  /  must  remain  a  normalized  distribution,  by 
definition.  Therefore,  probability  must  be  added  into 
the  volume  without  bias  towards  any  droplet 
characteristic  (in  proportion  to/). 

The  last  term  on  the  right  accounts  for 
redistribution  of  probability  in  a' -space  due  to  a 
gradient  in  the  mean  number  density  of  drops.  Suppose 
there  is  a  positive  gradient  in  the  mean  number  density 
of  drops  in  a  volume  and  the  drops  have  various 
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velocities.  Drops  that  are  moving  faster  than  the  mean 
velocity  will  leave  the  volume  in  greater  numbers  than 
those  that  are  moving  slower  than  the  mean.  Therefore, 
the  probability  of  finding  a  fast  drop  (faster  than  the 
mean)  decreases  while  the  probability  of  finding  a  slow 
drop  (slower  than  the  mean)  increases  and  thus  we  have 
a  shift  of  probability  in  a'  -space. 

We  now  have  two  equations  which,  together,  aie 
equivalent  to  equation  (3).  The  equation  for  the 
evolution  of  the  mean  number  density,  equation  (6),  is 
relatively  easy  to  solve  numerically  because  it  has  the 
form  of  the  wave  equation  for  which  there  are  many 
solution  techniques.  Conversely,  equation  (8)  presents 
us  with  many  of  the  same  problems  as  did  (3).  In  fact, 
equation  (8)  appears  more  complicated  than  (3).  Not 
only  do  we  still  have  to  solve  the  equation  across  all  of 
a -space,  but  we  have  introduced  two  additional  terms. 
However,  as  we  discuss  in  the  next  section,  we  aren’t  so 
much  interested  in  the  PDF  /  itself,  but  rather  various 
moments  of  that  function.  This  fact  will  be  used  to 
reduce  the  dimensionality  of  equation  (8). 

Moment  Transport  Equations 

One  of  the  primary  motivations  behind  this  study 
is  to  compute  various  statistical  averages  of  interest 
directly,  rather  than  post-average  volumes  of  simulation 
data.  In  this  section,  we  derive  a  set  of  moment 
transport  equations.  By  solving  these  equations,  we  are 
able  to  calculate  the  statistical  averages  at  any  point  in 
space  and  time  without  performing  a  numerical 
simulation.  Because  of  space  constraints,  we  only  give 
a  cursory  derivation  here.  A  detailed  derivation  can  be 
found  in  reference  [15]. 

Diameter  Moment  Transport  Equation 

Recall  that  the  expected  value,  or  ensemble  average, 
of  a  quantity  g(k)  is  defined  as 

{g)  =  j  g{k)f(k)dk  (9) 

k 

where  f(k)  is  the  PDF,  and  the  integral  is  over  the 
sample  space  of  k.  If  g(k )  is  a  positive,  integer-order 
power  of  k ,  then  (9)  defines  a  moment  of  f(k ) .  In 
particular,  the  average  drop  diameter  at  some  location  x 
and  t  is  defined  as 

(<p)  =  jj<t>f(<p,v;x,t)d\d<t)  (10) 

<t>  v 

In  the  same  manner,  we  can  take  moments  of  the 
transport  equation  for/.  Multiplying  equation  (8)  by  (j> 
and  recognizing  that  <p  is  an  independent  coordinate  so 
that  it  commutes  with  the  derivative  operators  gives 


^+v,(rv)+v,.(p)=  (n) 

0/  [  V,  •  (v) + (( v)  -  v)  •  V j  In  A] 

Integrating  over  all  diameter  and  velocity  space  yields 

M+Vj.(0v)  =  (0)Vj(v)  (i2) 

({0)(v)-(0v))-V,lnA 

This  equation  describes  the  evolution  of  the  mean  drop 
diameter  through  space  and  time.  By  making  use  of 
equation  (6)  and  the  product  rule,  (12)  can  be  cast  in 
conservative  form  as 

^M+Vi.(A(0v))  =  O  (13) 

Though  (13)  is  an  equation  for  the  transport  of  the  mean 
number  density  times  the  mean  drop  diameter,  rather 
than  just  the  mean  drop  diameter,  it  is  in  a  form  that  is 
more  convenient  than  equation  (12)  because  it  has  a 
form  similar  to  that  of  the  wave  equation.  Using  the 
same  procedure,  we  can  derive  transport  equations  for 
the  second  and  third  moments  about  the  origin  for  the 
drop  diameter: 

^l+Vx-(  A(*2v))  =  0  (14) 

&l  +  Vx-(^v))  =  0  (15) 

Velocity  Moment  Transport  Equations 

Before  we  can  derive  equations  for  the  drop  velocity 
moments,  we  must  make  some  changes  to  equation  (8). 
It  is  well  understood  that  droplet  velocities  correlate 
with  the  gas  velocities.  However,  /  does  not  contain 
any  information  about  the  gas.  Therefore,  we  must 
model  the  effect  the  gas  phase  has  on  the  drops.  To  do 
this,  we  divide  the  diameter  axis  into  equal-sized  “bins” 
and  derive  a  set  of  velocity  moment  equations  for  each 
bin.  We  can  then  consider  the  spray  as  the 
superposition  of  a  set  of  monodisperse  sprays.  This 
approach  is  similar  to  that  taken  by  Tambour11  when  he 
derived  a  series  of  sectional  equations,  each  section 
representing  one  bin.  Because  the  drops  do  not  vaporize 
or  breakup,  we  can  justify  treating  the  spray  in  this 
manner. 


4 

American  Institute  of  Aeronautics  and  Astronautics 


Using  conditional  probabilities,  we  split  /  into  two 
separate  probability  distributions  functions 

v ;  3c,  t)  =  p((j) ;  x ,  t)  q{ v ;  <j>,  x,  t)  (16) 

where  p((jr,x,t)d<j)  is  the  probability  that  a  drop  will  be 
within  a  diameter  interval  d<p  centered  about  0, 
conditioned  on  a  drop  existing  at  x  and  t ,  and 
q{v;<l>,x,t)dv  is  the  probability  that  a  drop  will  be 
within  a  velocity  element  dw  centered  about  v, 
conditioned  on  a  drop  with  diameter  <j>  at  x  and  t. 
Upon  substitution  of  (16)  into  (8),  we  have 

^+Mw*)+v.-(wO-  (17) 

m[V,(v)+((v)-v)-V>A] 


Now  let’s  focus  on  the  nth  bin  of  the  diameter  axis. 
Because  we  have  discretized  the  diameter  axis,  p($\x,t) 
must  also  be  discretized  and  treated  as  a  constant  in  each 
bin,  as  must  be  q(v\<j),x,t).  We  define  the  binned 
probability  density  f  ={pq)  as  the  probability  in  the 
bin  divided  by  the  width  of  the  bin 


<Pn 


(18) 


after  having  made  use  of  equation  (6). 

The  only  issue  that  remains  to  be  addressed  is  how 
to  handle  the  momentum  exchange  terms  on  the  right 
side  of  (21).  Recall  that  we  have  made  no  assumptions 
about  the  spray  that  limit  what  we  might  use  for  a  drag 
law,  allowing  us  to  specify  any  drag  law  of  interest 
without  significantly  changing  the  equations.  For  the 
moment,  we  will  assume  that  the  drops  decelerate  (or 
accelerate)  according  to  Stokes’  drag  law 


a  = — 


3ttu  0(v-w) 


n 


P ^ 


PtP 


(22) 


where  u  is  the  gas  velocity,  pg  is  the  gas  dynamic 
viscosity,  pg  and  p,  are  the  gas  and  liquid  densities, 
respectively,  and  vg-pg/p8  is  the  gas  kinematic 
viscosity.  Using  Stokes’  drag  is  not  a  limitation  of  the 
model.  Instead,  we  are  using  it  to  keep  the  analysis 
simple  while  we  examine  the  merit  of  the  overall  spray 
model.  Substituting  equation  (22)  into  (21)  gives 


=“((v*)„-k» 

Ln 


where 


We  also  define 


P  = 


jpd<j>  \pd<p 

± _ =  K _ 

J  d<p 


(19) 


We  can  thus  divide  equation  (17)  by  A  <pn  and  then 
integrate  over  the  width  of  the  diameter  bin  to  obtain 


~  +  ^ V,  •  {p'q  v)  ^ +  Vv  •  (p’qa)  = 
P  V  [  V,  •  (v) + ((v)  -  v)  •  V,  In  A] 


(20) 


Multiplying  (20)  by  the  drop  axial  velocity  v,  and 
integrating  over  all  velocity  space  in  the  nth  diameter 
bin  and  defining  Xn  =  Xp*  A(f)n  as  the  mean  number 
density  of  drops  in  the  ntb  diameter  bin,  we  arrive  at 


3(^(v^)J  +  v  (A  (-  )  \  =  2^{a  } 
dt  '  '  ,n)  A<j)n  '  ,n 


(21) 


pd)1 

_  r  IT  geo, n 

~  18Psys 


(24) 


is  the  representative  droplet  relaxation  time  and  <j)geo  n  is 
the  geometric  mean  diameter  of  the  bin. 

Similarly,  we  can  write  equations  for  the  mean- 
squared  axial  velocity,  mean  transverse  velocity,  mean- 
squared  transverse  velocity,  and  the  axial/transverse 
velocity  cross  moment 


(25) 


(26) 
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k 


(27) 


a(A»(V^)J  ,  y  ( 

l  (vv  V  \  ) 

dt  1  v 

n\  x  yin) 

-~L(2(V->V)’>„  “<va}„  ~(va)„) 


(28) 


Equations  (6),  (13)  through  (15),  (23),  and  (25)  through 
(28)  are  the  set  of  equations  that  we  will  use  to  describe 
the  evolution  of  the  spray  droplet  PDF.  However, 
because  we  can  treat  each  bin  as  a  separate  spray, 
equation  (6)  holds  in  each  bin,  using  each  bin’s  mean 
drop  velocity  (v)n  and  number  density  Xn .  Because  the 
overall  number  density  is  the  sum  of  the  binned  number 
densities  and  Xn=Xp*A<pn ,  we  can  reconstruct  the 
diameter  PDF  from  the  binned  number  densities.  This 
alleviates  our  need  to  solve  the  diameter  moment 
transport  equations  (13),  (14),  and  (15). 


Maximum  Entropy  Moment  Closure  (MEMO 

This  moment  method  is  based  on  a  hierarchical  moment 
structure  similar  to  that  encountered  in  turbulence 
modeling.  As  each  successive  moment  transport 
equation  is  derived,  at  least  one  new,  higher-order 
moment  is  introduced.  These  unclosed  moments  can  be 
evaluated  by  making  use  of  the  Maximum  Entropy 
Formalism  (MEF).  In  his  study  of  communications 
theory,  Shannon18  developed  the  concept  of  information 
entropy  as  a  measure  of  probabilistic  uncertainty.  For  a 
continuous  distribution,  this  measure  is  given  by 

S  =  -£  f(x)  In  f{x)dx  (29) 


Later,  Jaynes19  showed  that  an  infinite  number  of  PDFs 
are  consistent  with  a  set  of  known  constraints,  but  that 
the  one  that  should  be  chosen  is  the  one  with  maximum 
entropy.  If  a  PDF  with  less  entropy  (i.e.,  less 
uncertainty)  were  used,  it  would  imply  the  existence  of 
some  additional  knowledge.  However,  since  all  the 
available  knowledge  was  applied  in  the  form  of 
constraints,  no  additional  knowledge  would  exist.  Thus 
it  would  be  inappropriate  to  choose  any  PDF  other  than 
the  one  with  maximum  entropy.  This  PDF  is  the  most 
unbiased  distribution  possible  within  the  given 
constraints. 

Because  they  have  the  greatest  influence  on  the 
shape  of  the  PDF,  we  will  constrain  the  lowest-order 


moments  (i.e.  means,  variances,  etc.).  These  moments 
correspond  to  those  for  which  we  have  derived  transport 
equations.  Equation  (29)  is  maximized  by  way  of  the 
method  of  Lagrange  multipliers.  This  leads  to  a  PDF 
of  the  form 


/-exp 


-Po-ll 


(30) 


where  the  /J’s  are  the  Lagrange  multipliers,  gr(a ')  is 
the  r*h  functional  quantity  to  be  averaged,  and  M  is  the 
number  of  moment  constraints.  The  coefficients  are 
determined  from  a  non-linear  system  of  coupled 
differential  equations  which,  in  general,  must  be  solved 
numerically. 

Sellens20  and  Chin  et.  al.21  have  made  use  of  the 
MEF  in  their  work  on  predicting  droplet  distributions 
resulting  from  breakup  of  a  liquid  sheet  or  jet.  In  this 
work,  the  breakup  process  is  considered  independently 
from  gas-phase  constraints,  and  it  is  not  their  intention 
to  solve  for  the  evolution  of  the  droplet  distribution 
throughout  space  and  time,  but  rather  to  obtain  a 
distribution  for  all  the  drops  in  the  spray.  In  contrast, 
our  work  is  focused  on  describing  the  evolution  of  a 
droplet  PDF  at  every  point  in  the  flow  domain  and  how 
that  evolution  is  influenced  by  and  interacts  with  the 
gas  phase.  In  our  approach,  we  model  the  evolution  of 
the  complete,  coupled  spray  flow  using  moment 
equations.  The  MEF  is  introduced  as  a  method  to  close 
those  equations. 

Quasi-One-Dimensional  Spray  Problem 

To  explore  the  usefulness  of  this  approach,  we 
make  use  of  a  quasi-one-dimensional  spray  flow.  The 
purpose  of  this  is  to  pose  a  spray  problem  that  is 
understandable,  physically  interesting,  and  tractable. 
Moreover,  we  wish  to  look  at  a  simple  problem  to 
determine  if  there  are  any  flaws  in  the  proposed  closure 
models  which  will  present  themselves  more  readily  in 
this  problem  than  in  a  more  complicated,  multi¬ 
dimensional  problem. 

Figure  (1)  shows  the  geometry  for  this  problem. 
At  x  =  0  there  is  a  two-dimensional  array  of  spray 
injectors  that  extends  infinitely  in  the  y-direction.  Each 
injector  has  been  placed  at  random  on  this  array  and 
points  in  the  jc-direction,  delivering  a  distribution  of 
drops  into  an  incompressible  gas,  which  initially  has  a 
low  level  of  turbulence.  The  array  and  each  injector  on 
it  extend  infinitely  in  directi  on,  resulting  in  a  series  of 
wedge  sprays  of  spherical  drops. 

Now  imagine  an  ensemble  of  these  arrays  where, 
on  each,  the  injectors  have  been  placed  at  random. 
Because  we  wish  to  develop  a  quasi-one-dimensional 
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flow,  we  restrict  the  PDF  for  the  drop  velocities  at  each 
injector  to  have  a  mean  transverse  (^-direction)  velocity 
of  zero.  If  there  is  an  infinite  number  of  realizations  in 
the  ensemble,  then  there  is  no  preferential  location 
along  the  transverse  axis  since  there  is  equal  probability 
for  a  drop  to  have  a  positive  y-velocity  of  some 
magnitude  as  there  is  for  a  drop  to  have  a  negative  y- 
velocity  of  the  same  magnitude.  Similarly,  we  have  a 
zero-mean  gas  velocity  component  in  the  y-direction, 
though  fluctuations  do  occur.  Because  there  is  no 
preferential  location  along  the  transverse  axis,  there  arc 
no  variations  along  that  axis  across  the  ensemble ,  and 
therefore,  derivatives  of  averaged  quantities  in  the  y- 
direction  are  zero.  So,  the  three  primary  constraints  we 
have  for  this  problem  are 

(“,)  =  «■  (v,)  =  0,  (31) 

Also,  so  that  no  mean  flow  develops  in  the  transverse 
direction,  there  can  be  no  correlation  between  the  x-  and 
y-components  of  velocity  in  either  phase.  Thus,  all 
cross-component  and  cross-component/cross-phase 
moments  (i.e.  {yxvy)’{uxvy)’  etc-)  eclual  zera 
Substituting  these  constraints  into  the  velocity  moment 
transport  equations  yields 


d{Uj)_  d(p) 
dt  dx 


^fl+3nvt^yx-ux))  (37) 


Equation  (36)  says  that  the  mean  gas  velocity  is  a 
constant  in  space.  If  we  specify  the  inlet  mean  gas 
velocity  as  a  constant  in  time,  then  the  mean  gas 
velocity  throughout  the  domain  is  constant  in  both 
space  and  time.  This  makes  solving  the  gas  momentum 
equation  unnecessary  unless  one  is  interested  in  the 
pressure.  Since  our  present  interests  are  in  the  discrete 
phase  and  not  the  gas,  we  do  not  need  to  solve  equations 
(36)  and  (37).  However,  a  k-e  turbulence  model, 
modified  to  account  for  the  drops,3  is  solved  since  the 
gas-phase  turbulence  is  incorporated  into  the  cross-phase 
moments  appearing  in  equations  (34)  and  (35). 

The  moment  equations  are  solved  using  an 
implicit,  first-order,  upwind  scheme  for  the  spatial 
derivatives  and  a  second-order  backwards  difference  for 
the  time  derivatives.  A  first-order  scheme  was  chosen 
for  the  spatial  derivatives  due  to  the  difficulties 
associated  with  solving  moment  equations.  Moment 
variables  have  relations  among  them  that  must  be 
preserved,  constraining  their  values.  For  example,  the 
relation  between  the  first  and  second  moment  of  a  PDF 
is 


(x2)-(x)2>0  (38) 


9A„  [  d(*.(v*),)_0 


dx 


(32) 


(33, 


dt 

Ln 

3K(V*).)  ,  d(A«(V*V')„) 


dt 


dx 


(34) 


(35) 


If  constraints  (31)  are  applied  to  the  incompressible 
RANS  equations,  we  obtain 


3(»x) 

3* 


=  0 


(36) 


where  X  is  any  random  variable.  If  this  relation  is  not 
maintained,  then  the  numerical  solution  may  either 
become  unstable  or  give  non-physical  answers.  Second- 
order  numerical  solutions  of  hyperbolic  equations  can 
cause  oscillations  in  the  solution  in  regions  of  strong 
gradients,  even  when  artificial  dissipation  is  added. 
Being  non-physical,  those  oscillations  could  cause  the 
numerical  solution  of  the  variance  to  become  negative, 
corrupting  the  remainder  of  the  solution. 

For  comparison,  a  Lagrangian  simulation  of  the 
quasi-one-dimensional  spray  is  performed.  Individual 
drops  are  tracked  from  the  injector  along  their 
trajectories.  Because  the  mean  gas  velocity  is  a 
constant  in  space  and  time,  the  gas  momentum  and 
continuity  equations  are  not  solved,  though  a  k-e 
model  is  used  to  simulate  the  turbulence,  giving  each 
drop  a  random  “kick”  over  a  duration  equal  to  the  eddy 
turnover  time  or  the  drop  residence  time,  whichever  is 
shorter. 


Results 

Table  1  lists  the  moments  used  to  specify  the 
droplet  PDF  at  the  injector  and  the  gas  phase 
characteristics.  The  diameter  PDF  is  shown  in  Figure  2 
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and  the  velocity  PDFs  are  Gaussian.  The  inlet 
condition  for  the  turbulent  kinetic  energy  corresponds  to 
a  RMS  value  of  five  percent  of  the  mean  gas  velocity 
and  the  rate  of  dissipation  of  turbulent  kinetic  energy  is 
£  =  /:/100sec. 

Figure  3  shows  profiles  of  the  mean  number 
density.  It  is  clear  where  the  leading  edge  of  the  spray 
is  at  each  time  interval  and  how  it  propagates  through 
the  domain  as  a  concentration  wave.  The  mean  number 
densities  behind  it  are  much  greater  than  the  injection 
mean  number  density.  This  is  what  we  would  expect 
since  the  drops  are  injected  into  a  slower  gas.  The  drops 
rapidly  decelerate  and  accumulate  just  behind  the  leading 
edge.  This  accumulation  of  drops  continues  until  the 
drops  have  been  slowed  to  the  mean  gas  velocity.  The 
results  compare  well  with  the  Lagrangian  simulation  in 
the  region  behind  the  leading  edge.  At  the  leading  edge, 
the  MEMC  results  deviate  from  the  simulation. 
However,  this  is  due  to  the  numerical  dissipation 
associated  with  the  first-order  numerical  scheme. 

At  the  leading  edge,  the  mean  number  density  is 
quite  low.  When  drops  first  arrive  at  some  spatial 
location,  only  a  few  are  present,  those  that  were  large 
enough  not  to  have  been  greatly  influenced  by  drag. 
This  is  shown  in  Figure  4  where  the  expected  drop 
diameter  is  plotted.  At  the  leading  edge,  there  is  a 
greater  mean  drop  diameter.  This  is  caused  by  the 
strong  slip  velocity  present  in  the  flow,  decelerating  all 
but  the  most  massive  drops.  As  we  move  upstream,  we 
encounter  the  smaller  drops  that  were  quickly 
decelerated. 

Figure  5  shows  how  the  standard  deviation  of  the 
diameter  PDF  evolves.  At  the  leading  edge,  the  standard 
deviation  is  large.  This  comes  from  the  fact  that  there 
are  so  few  drops  out  on  the  wave  front.  Those  drops 
present  are  large  due  to  the  developed  size-velocity 
correlation.  Figure  2  shows  that  there  is  a  wide  range 
of  large  diameters  that  those  drops  have.  It  can  be 
argued  that  the  shape  of  the  standard  deviation  curve 
behind  the  leading  edge  is  a  phenomenon  that  arises 
from  the  spatial  spreading  of  the  drops  that  were 
initially  injected,  however,  because  the  numerical 
diffusion  at  the  leading  edge,  this  effect  may  appear 
more  significant  than  it  really  is. 

Figure  6  shows  the  mean  drop  axial  velocity 
profiles.  The  region  closest  to  the  injector  is  where  the 
highest  velocity  slip  occurs.  There  are  a  number  of 
drops  passing  through  this  region  that  are  moving  faster 
than  the  mean  gas  velocity.  Due  to  their  lower  masses, 
most  of  those  drops  rapidly  decelerate,  reducing  the 
mean  velocity.  The  mean  velocity  of  these  drops 
approaches  the  mean  gas  velocity  within  3  meters.  At 
the  leading  edge,  however,  the  mean  velocity  again 


begins  to  increase.  To  have  reached  the  leading  edge, 
those  drops  had  to  be  the  fastest  coming  out  of  the 
injector,  and  because  they  are  also  some  of  the  largest, 
they  haven’t  been  affected  by  drag  as  much  as  the 
smaller  drops.  Thus,  the  mean  velocity  is  higher  at  the 
edge.  As  the  wave  front  continues  to  move  forward, 
those  large  drops  slow  as  evidenced  by  the  decrease  in 
the  mean  velocity  on  the  front  as  time  goes  on. 

Figures  7  and  8  show  the  standard  deviations  of  the 
axial  and  transverse  velocity  distribution  functions. 
They  are  measures  of  the  width  of  the  velocity  PDFs. 
The  larger  widths  occur  near  the  injector  because  the 
drops  have  not  been  influenced  by  the  gas  phase  long 
enough  to  have  approached  the  gas  velocity. 
Downstream,  the  standard  deviations  rapidly  decay  until 
we  reach  the  wave  front.  It  will  take  a  longer  time  for 
the  larger  drops  to  reach  the  gas  velocity,  but  as  the 
leading  edge  propagates,  the  standard  deviations  become 
smaller  and  smaller.  Again,  there  is  good  agreement 
with  the  simulation. 

One  thing  of  note  is  the  magnitude  of  the  drop 
axial  velocity  standard  deviation  relative  to  the  mean 
drop  velocity.  At  the  injector,  the  standard  deviation 
was  prescribed  as  one-tenth  of  the  injection  mean.  In 
Figure  7,  however,  we  see  it  is  much  wider  in  the 
region  of  strong  slip.  This  results  from  the  fast  drops 
quickly  moving  ahead,  and  slower  ones  lagging.  Many 
drops  close  to  the  injector  have  decelerated,  especially 
the  smallest  ones.  There  is  also  a  significant  number  of 
larger  drops  in  this  region  that  are  moving  with 
velocities  closer  to  the  mean  injection  velocity  or  faster. 
Most  of  these  drops  are  the  more  massive  ones  not 
having  been  greatly  influenced  by  drag  as  yet.  Thus,  in 
the  region  of  strong  slip,  there  is  a  much  wider 
distribution  of  velocities  as  a  result  of  deceleration. 
This  is  the  reason  for  the  much  higher  axial  velocity 
standard  deviation. 

The  simulation  results  in  the  figures  are  the 
averages  of  ten  runs.  Because  of  the  nature  of  the 
simulation,  and  because  the  results  are  averages  of  only 
a  few  realizations,  there  is  a  great  deal  of  noise  in  the 
simulation  data.  To  smooth  this  out,  and  to  better 
approximate  an  ensemble  average,  we  would  ideally  like 
to  average  over  thousands  of  simulations.  However, 
because  of  the  computational  time  required,  we  are 
prohibited  from  doing  this.  The  best  we  can  do  is  time 
average,  but  we  are  limited  to  doing  this  only  in  regions 
that  have  reached  a  steady  state.  We  are  not  able  to 
capture  the  averaged  dynamics  of  the  leading  edge  using 
a  simulation  unless  we  ensemble  average.  This  is  a 
significant  advantage  that  the  current  model  has  over  the 
Lagrangian  simulation. 

A  second  advantage  is  the  savings  in  computational 
time  for  getting  detailed  statistics.  Because  we  always 
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have  an  approximation  to  the  PDF  as  part  of  the 
solution  to  the  MEMC  model,  we  can,  at  any  time, 
integrate  that  PDF  and  obtain  any  statistics  of  interest. 
This  is  not  true  with  the  Lagrangian  simulation,  since, 
as  the  level  of  statistical  detail  increases,  so  does  the 
time  required  to  collect  sufficient  data.  Figure  9 
compares  the  nondimensional  times  required  to  obtain 
different  types  of  statistics  between  the  two  models  for 
this  problem.  Other  problems  may  have  different  timt 
results,  though  the  general  trends  will  be  the  same.  The 
MEMC  calculation  was  run  for  six  seconds  of  flow 
time.  Bar  A  represents  the  time  required  to  run  the 
simulation  for  the  same  amount  of  flow  time  and  obtain 
a  “snapshot”  of  the  simulation  at  that  point.  Bar  B 
denotes  the  amount  of  time  required  to  get  steady-state 
statistics  at  all  grid  locations  up  to  100  mm.  Bars  C 
and  D  show  the  time  needed  to  obtain  more  detailed 
steady-state  statistics  up  to  100  mm  downstream  of  the 
injector.  In  case  C,  statistics  were  computed  only  on 
drops  with  a  diameter  greater  than  100  (im  and  in  case 
D,  only  on  drops  with  a  diameter  greater  than  100  pm 
and  a  velocity  greater  than  23  mm/s.  Note  that  as  the 
level  of  detail  increases,  so  does  the  simulation  time. 
To  obtain  averaged  data  from  the  simulation  for  the 
transient  region  requires  an  even  greater  amount  of  time 
since  we  must  ensemble  average  this  region,  requiring 
many  simulation  runs. 

Conclusions 

We  have  developed  a  model  which  describes  the 
evolution  of  a  spray  flow  by  transporting  a  series  of 
moments  of  the  droplet  probability  density  function  and 
closing  those  equations  with  a  maximum  entropy 
model.  The  model  was  tested  on  a  quasi-one- 
dimensional  spray  and  results  for  several  quantities  of 
the  spray  were  obtained,  including  the  expected  number 
density,  expected  droplet  diameter,  and  expected  droplet 
velocity.  We  obtained  acceptable  agreement  with 
statistical  data  obtained  from  a  Lagrangian  simulation, 
both  in  the  steady-state  regions  and  the  transient  region. 
We  also  showed  that  the  present  model  can  offer  a 
significant  computational  time  saving  over  a  simulation 
when  detailed  statistics  are  desired. 
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Table  1  Injector  Conditions 


Figure  2  Injector  diameter  PDF 


Expected  Drop  Number  Density 


Figure  3  Mean  drop  number  density 
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Figure  4  Mean  drop  diameter 


Figure  7  Drop  axial  velocity  standard  deviation 


Drop  Diameter  Standard  Deviation 


Drop  Transverse  Velocity  Standard  Deviation 
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lure  5 


Drop  diameter  standard  deviation 


Figure  8  Drop  transverse  velocity  standard  deviation 


Expected  Drop  Axial  Velocity  Statistics  Collection  Times 


Position  (mm) 

Figure  6  Mean  drop  axial  velocity  Figure  9  Simulation  vs.  MEMC  computational 

times. 
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